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ABSTRACT 

Direct  integration  leads  to  the  potential  gradient  of  a  rectangular  plate.  Lines  with 
equal  spacing  form  a  grid  for  integration  over  the  infinite  plane  A  source  distribution 
is  approximated  by  pattern  functions  which  are  centered  over  each  grid  point.  The 
pattern  functions  are  based  upon  the  sine  quotient  function.  Fourier  analysis  replaces 
integration  in  physical  space  with  integration  in  wave  number  space  The  range  of 
integration  is  limited  to  lines  of  finite  length  in  wave  number  space.  The  computation 
of  potential  gradient  is  provided  by  subroutines. 


n 


INTRODUCTION 

The  flow  around  a  ship  is  partly  a  radial  flow  from  a  source  distribution  and  partly 
a  circular  flow  around  a  vortex  distribution.  In  either  case  the  velocity  is  the  integral 
of  the  product  of  an  inverse  square  and  the  density  of  sources  or  vortices.  The  inverse 
square  contributes  a  spike  to  the  integrands  and  numerical  integration  is  not  feasible. 
A  technique  for  integration  through  the  spike  may  be  found  from  an  analysis  of  the 
analogous  problem  in  the  flow  over  a  plane. 

The  potential  and  the  potential  gradient  of  a  source  or  charge  which  is  distributed 
over  a  plane  have  applications  in  hydrodynamics,  electromagnetism,  and  gravitation. 
The  field  of  an  irregular  three-dimensional  source  distribution  can  be  constructed 
with  source  distributions  over  each  of  a  stack  of  planes.  The  potential  and  the  potential 
gradient  at  a  field  point  near  a  plane  can  be  computed  by  two  methods. 

In  a  direct  method,  the  potential  at  a  field  point  is  the  integral  over  the  plane  of 
the  product  of  inverse  distance  and  the  source  density,  while  the  potential  gradient 
is  the  integral  over  the  plane  of  the  product  of  an  inverse  square  of  distance  and  the 
source  density.  Wherever  there  is  a  discontinuity  in  source  density  the  potential 
gradient  is  infinite 

In  a  Fourier  method,  the  potential  is  expressed  as  a  solution  of  Laplace's  equation, 
and  the  solution  is  adjusted  to  fit  the  boundary  conditions.  The  solution  is  the  sum 
of  products  of  complex  exponential  functions,  and  the  source  distribution  is 
approximated  as  a  Fourier  series.  A  Fourier  transform  replaces  integration  in  physical 
space  with  integration  in  wave  number  space  The  time  of  computation  depends  upon 
the  distance  to  the  field  point.  There  are  two  methods  for  the  integration  in  wave 
number  space.  In  the  first  method,  Gauss  integration  multipliers  are  applied  in  a 
roving  interval.  The  time  increases  with  distance.  In  the  second  method,  integration 
by  parts  leads  to  recurrence  equations.  The  time  is  independent  of  distance. 

DIRECT  INTEGRATION 

Let  x,y,z  be  Cartesian  coordinates  in  a  right-handed  coordinate  system  with  x  and  y 
in  the  plane  and  with  z  above  the  plane.  Let  o(x,  y)  be  the  source  density  at  a  point 
with  coordinates  x,  y.  Then  the  potential  is  given  by  the  equation 

.  ^       C   C  a(u- v)  ,    \ 

<p(x,y.z)  =  dudv  (1) 

J  J    V(x  -  u)2  +  (y  -  v)z  +  zz 
where  x,y,z  are  coordinates  of  a  field  point  and  u,  v  are  coordinates  of  a  source  point. 

RECTANGULAR  PLATE 

Let  a  rectangular  plate  be  placed  on  the  plane.  Let  the  length  of  the  plate  be  2a, 
and  let  the  breadth  of  the  plate  be  26.  Let  the  origin  be  at  the  center  of  the  plate 
with  the  x-axis  in  the  direction  of  length  and  with  the  y-axis  in  the  direction  of 
breadth. 

The  potential  of  a  plate  with  unit  source  density  is  given  by  the  equation 

.+b-V     pa-x  dudy 

<p(x,y,z)  =  ,   „        „        =  (2) 


/.-«,-„    pa-x  dud 

>(z.2/.  z)  =  /    - 

J  -b-y    J-a-x     VU     +  V 


Initial    integration    leads    to    inverse    hyperbolic    functions,    then    final    integration    is 
completed  with  integrations  by  parts. 


The  potential  is  given  by  the  equation 

</>  =  (a  -  x)sinh        ,  +  (a  -  x)sinh 


V(a-x)2  +  22  V(a  -  x)2  4-  z2 

(6  +  2/) 


+  (a  +  x)sinh   '     ,  +  (a  +  x)sinh   ! 

V(a  +  x)2  +  zz  v(a  +  x)2  +  z2 

.  _.         (a  -  x)  _.         (a  +  x) 

+  (b  -  y)sinh        ,  +  (b  -  y)sinh 


V(6-y)2+  z2  V(6-Ty)2  +  Z2 


+  (b  +  y)sinh   '     ,  +  (b  +  y)sinh 

V(6  +  y)2  +  22  n/(6T^)2  +  22 

(a-x)(b-y)  _  (a  -  x)(b  +  y) 

-  2  tan    ' -  z  tan 


2V(a-x)2+  (6-y)2  +  22  2V/(a-i)2  +  (6  +  2/)2  +  22 

(a  +  x)(6-y)  (a  +  x)(6  +  y) 

-  2  tan   ' — ,  = 2  tan     — -. — (3) 

2V(a  +  x)2  +  (b  -  y)2  +  z2  2V(a  4  x)2  +  (6  +  y)2  +  22 

Partial  differentiation  and  cancellation  of  terms  gives  the  components  of  the  gradient 
of  the  potential. 

Differentiation  with  respect  to  x  leads  to  the  equation 

dp        .    ,_,        (b-y)  ,        (b  +  y) 

=  sinh         .  +  sinh         ,  — 

dx  Vi^T-  x)2  +  22  V(a-x)2  +  22 

.    uM         (b-y)  .      _, (b  +  y)  m 

-sinh        ,  -  sinh        ,  (4) 

V  (a  +  x)2  +  22  V(a  +  x)2  +  22 

differentiation  with  respect  to  y  leads  to  the  equation 

dip        .         ,         (a  -  x)  (a  +  x) 
=  sinh        ,                         +  sinh     - 


dy  V(6-y)2  +  22  N/(6_:y)2  +  22 

u-i         (a_:r)  ■    ,-i         (a  +  x) 

-  sinh         ,  -  sinh    '     ,  (5) 

V(6  +  y)2+22  V(6  +  y)2  +  22 

and  differentiation  with  respect  to  z  leads  to  the  equation 

d<p  (a-x)(b-y)  _,  (a  -  x)(b  +  y) 

=  tan     — ,  +  tan   '       , 

dz  2V(a-  x)2+  (6  -  y)2  +  22  zV(a  -  x)2  +  (6  +  y)2  +  22 

,*       _,  (a  +  x)(b-y)  _  (a  +  x)(b  +  y) 

+  tan     — ,  +  tan   1 — ,  (6) 

z\/(a  +  x)2  +  (b  -  y)2  +  22  2V(a  +  x)2  +  (b  +  y)2  +  22 

These  equations  were  published  in  a  previous  report1.  They  are  used  in  the  following 
subroutines. 


SUBROUTINE   RPLTP   (AA,  AB,  AX,  AY,  AZ,   FP) 

********************************************************************************************* 

FORTRAN  SUBROUTINE  FOR  POTENTIAL  OF  RECTANGULAR  PLATE 

********************************************************************************************* 

The  half-length  a  of  the  plate  is  given  in  argument  AA,  the  half-breadth  b  of  the 
plate  is  given  in  argument  AB,  and  the  Cartesian  coordinates  x,  y,  z  of  the  field  point 
are  given  in  arguments  AX,  AY,  AZ.  The  potential  of  the  plate  is  stored  in  function  FP. 

SUBROUTINE   RPLTG  (AA,  AB,  AX,  AY,  AZ,   FX,   FY,   FZ) 
********************************************************************************************* 

FORTRAN  SUBROUTINE  FOR  POTENTIAL  GRADIENT  OF  RECTANGULAR  PLATE 

*************************************************************************************** 

The  half-length  a  of  the  plate  is  given  in  argument  AA,  the  half-breadth  6  of  the 
plate  is  given  in  argument  AB,  and  the  Cartesian  coordinates  x,  y,  z  of  the  field  point 
are  given  in  arguments  AX,  AY,  AZ.  The  components  of  the  gradient  are  stored  in  functions 

FX,  FY,  FZ. 

NONUNIFORM  PLATE 

A  major  achievement  was  the  integration  when  the  source  density  was  expressed 
as  power  polynomials  in  the  surface  coordinates.  The  computation  of  potential  required 
a  triplex  of  recurrence  relations.  The  computation  is  provided  by  the  subroutine  RPLTM. 
The  analysis  has  been  presented  in  a  previous  report1.  That  RPLTG  and  RPLTM  give  the 
same  gradients  has  been  verified  by  computation. 

FOURIER  ANALYSIS 
The  Fourier  transform  for  a  function  on  a  plane  is  given  by  the  equations 


A(a,p)  =  -^\    \  F(x,  y)  e-i{ax+Py)  dx  dy 
F(x,y)  =   I    J  A{a,p)ei(ax+Iiy)  da 


(7) 


d/8  (8) 

where  x,  y  are  Cartesian  coordinates  in  physical  space,  a,/?  are  Cartesian  coordinates 
in  wave  number  space,  F(x,  y)  is  a  function  in  physical  space,  and  A(a,  /?)  is  the  Fourier 
amplitude  in  wave  number  space.  A  potential  (p(x,y,  z)  of  a  source  distribution  on  the 
plane  can  be  constructed  with  the  equation 


<p(x,y,z)=   j    (A(a,(3)e-^F^M+i{ax+Pv)da 


d/S  (9) 


which  gives  a  solution  of  Laplace's  equation  wherever  2*0.  The  vertical  component 
of  the  gradient  at  z  =  0  is  given  by  the  equation 


~dz 


Va2  +  /?2yl(a,  0)  eiiax+fiv)  da  dp  (z  -*  0)  (10) 


where  the  sign  is  +  if  z  >  0  and  the  sign  is  -  if  z  <  0.  In  accordance  with  the  Gauss 
theorem  the  difference  in  -  dcp/dz  on  opposite  sides  of  the  plane  is  4na(x,  y)  where 
a(x,  y)  is  the  surface  source  density  at  the  plane.  Thus  the  Fourier  amplitude  for  an 


arbitrary  distribution  of  source  density  is  given  by  the  equation 

A(a,fi)  =  }         ?    f  \  a{x,y)e-i(ax^y)dxdy  (11) 

The  range  of  integration  of  density  is  the  infinite  plane. 

RECTANGULAR  PLATE 

Let  a  rectangular  plate  be  placed  on  the  plane.  Let  the  length  of  the  plate  be  2a, 
and  let  the  breadth  of  the  plate  be  26.  Let  the  origin  be  at  the  center  of  the  plate 
with  the  x-  axis  in  the  direction  of  length  and  with  the  y-axis  in  the  direction  of 
breadth.  Let  u,  v  be  the  coordinates  of  a  source  point  on  the  plate  and  let  x,  y,  z  be 
the  coordinates  of  a  field  point  above  the  plate. 

The  potential  of  a  plate  with  unit  source  density  is  given  by  the  equation 


<fi(x.y,  z)  = 


i,un 


■N/a2+/J2|z|  +  i|a(i-u)+/?(y-v)| 


dadpdudv  (12) 


Va2  +  /?2 
Completion  of  the  integration  leads  to  the  equation 

,       2    f   f  sin(aa)  sin(60)  e-\fi^il*l+*(«*+*v) 

<p(x.y.z)  =  -\        > - ^  .  da  dp  (13) 

7T  J   J  a  0  Va2  +  02 

What  is  interesting  about  this  equation  is  that  the  distribution  of  source  density  over 
a  finite  plate  introduces  the  factors  sin(aa)/a  and  sin(6/!?)//3  into  the  Fourier  amplitude. 
The  factors  concentrate  the  spectrum  of  the  Fourier  integral  into  bands  which  are 
centered  with  respect  to  the  line  spectrum  of  a  trigonometric  polynomial. 

INTERPOLATION 

The   coefficient   function    for   Lagrange   interpolation   in    an    infinite   set   of   equally 
spaced  data  is  just  the  infinite  product  for  the  sine  quotient  function  in  the  equation 


sin  u 


u 


■JJ.('-H?)  (14) 


It  is  equal  to  unity  at  the  origin  and  is  equal  to  zero  wherever  else  u  is  a  multiple 
of  7T.  It  is  analytic  and  diminishes  with  distance  from  the  origin. 

An  arbitrary  function  f(u)  is  expressed  in  terms  of  coefficient  functions  with  centers 
at  kn  by  the  equation 

.    ,      kZX°°    .      ,  sin(-u  -  kn) 

f(u)=    E    f(uk) K—- '-  (15) 

t=-oo  u-k-n 

where  f{uk)  is  the  discrete  value  of  f(u)  at  the  fcth  node. 

The  product  of  functions  which  are  centered  at  kn  and  mn  is  given  by  the  function 

sin(u  -  kn)  sin(u  -  mn) 

i     .  v — r v~  (16) 

(U  -  klT)  [U  -  7717T) 

Resolution    into    partial    fractions    and    application    of    the    addition    theorem    for 


trigonometric  functions  leads  to  the  equation 

f+0°  sm(u-kTr)  sm(u-mn)  (-  l)*+TrT   f+°°  sin2(u-k-n)  j  f+0°  sin2(u-m7r)  j     1,      . 

— — — -  du  =  — —  ■ du  -  dit    (17) 

J.^      (u-A;7t)         (u-mn)  (fc-m)7rL  J_oo         u-kn  J _„         u-mn 

Inasmuch  as  the  integrands  are  odd  their  integrals  are  zero  and  the  coefficient 
functions  are  orthogonal  with  respect  to  integration  They  may  be  used  as  the  basis 
for  approximations  by  least  squares. 

An  important  property  is  expressed  by  the  equation 


r 


sin  au 

du  =  <, 

u 


—  7T  if  a  <  0 
Oifa  =  0  (18) 

I  +n  if  a  >  0 


where  a  is  a  constant  of  arbitrary  magnitude   Integration  by  parts  leads  to  the  equation 

J  +  oo               2                              n  +  oo                0 
sin'u             I        sin  2u 
T-du=  du  =  7T  (19) 
-co       u*                  J_x         u 

which  gives  the  normalization  factor  for  least  squares  approximation. 

SOURCE  PATTERN 
A  universal  source  pattern  for  the  infinite  plane  is  given  by  the  equation 


.    (txx\    .    (ny\ 
sin   —    sin   — 
\2o/        \26/ 


F{x,y)  =  7 —. \ (20) 

■nx\(  rty 


2a/\2b 

The  function  F(x,y)  is  unity  at  the  origin  and  is  zero  on  all  equally  spaced  nodal  lines 
which  do  not  pass  through  the  origin  The  spacing  between  nodal  lines  is  2a  in  the 
x-direction  and  is  2b  in  the  z/-direction.  The  function  is  analytic  and  diminishes  with 
distance  from  the  origin.  Application  of  the  Fourier  transform  and  use  of  the  Euler 
theorem  introduces  the  difference  of  two  cosines  and  introduces  the  sum  of  two  sines 
into  the  integrand  of  the  Fourier  amplitude.  The  difference  between  the  cosines  is 
zero  at  the  origin  and  is  an  even  function  otherwise.  The  difference  between  the  cosines 
is  eliminated  by  integration.  There  remains  the  sum  of  sines  as  expressed  by  the 
equation 

w      ^  1       f  f  sin(£-a):r  +  sin(^  +  a)x  sin(£-j3)y  +  sin(^+j8)y 

4(a,0)  =  —- g dxdy  (21) 

167T3  J   J  /7rx\  /ny\ 

\2a/  \26/ 

The  sines  cancel  each  other  during  integration  wherever  a,  (3  are  outside  a  rectangle 
in  wave  number  space.  The  interior  of  the  rectangle  is  located  where  a,/S  satisfy  the 
inequalities 

7T  7T  7T  7T 

^a^  +  —  gfig  +  —  (22) 

2a  2a  26  26 


Thus  the  function  F(x,  y)  is  given  by  the  equation 


a6     P    2b     P     2a       ..  „    . 

F(x,y)  =  —^\  ex{ax+/3y)  da  d/3  =  a(x,y)  (23) 

~2b         ~2a 


and  the  potential  (p(x,  y,  z)  is  given  by  the  equation 


2ab    C  +  ib   C  +  t  e-N/a2+/>2l*!+i<a*^y> 
<p(x,y,  z)  =  ~j    n    J    n  ,    „  .    „, dadP 


2b  2a 


Va2  +  0: 


Further  integration  is  possible  after  Cartesian  coordinates  a,/?  are  replaced  by  polar 
coordinates  /c,  6.  The  potential  is  expressed  by  the  equation 

<p(x,y,  z)  =  —    I     (  e_,c|*l+i*(a:cose+v,in')d/cd0  (25) 

Then  integration  with  respect  to  /c  replaces  surface  integration  throughout  the  interior 
of  the  rectangle  with  line  integration  along  the  perimeter  of  the  rectangle. 

The  three  components  of  the  gradient  are  given  by  the  same  formulations  except 
for  a  weighting  function  w.  Thus  for  the  components 


(26) 


d<p                                                      dip  dtp 

dx                                                    dy  dz 

the  values  of  w  arc 

-  rcosfl  -  i  sin  6  +1                   (27) 


Then  each  component  is  the  sum 


A+B  (28) 


of  two  integrals  where  A  is  the  integral  along  vertical  sides  of  the  rectangle  and  B  is 
the  integral  along  horizontal  sides  of  the  rectangle. 

For  the  vertical  sides  the  variables  are  related  in  accordance  with  the  equations 

*  =  tan0  (29) 

t                                                        1                                                      dt  ,      ' 

sin  6  =  —=  cosg  =  —.  dd  =  +  =         (30) 

7T 

a  =  k  cos  9  -  —  (31) 

2a 

Let  a  parameter  <5  be  defined  by  the  equation 

6  =  £-  i|2|Vl  +  t*  -  i{x  +  yt)\  (32) 

2a 


Then  the  integral  A  is  given  by  the  equation 


6    f+b  1  -  (1  +  6)e~6  .      , 

A  =  tt-\ — ~ wdt  (33) 

b 


for  integration  along  ^oth  vertical  sides. 


For  the  horizontal  sides  the  variables  are  related  in  accordance  with  the  equations 

*  =  cot0  (34) 

1                                                         t                                                      dt  ,      . 

sin0=^=  cos0=^==  d6  = -^         (35) 

(36) 
Let  the  parameter  6  be  defined  by  the  equation 

6  =  £-  \\z\y/l  +tz  -  i(xt  +  y)\  (37) 

CO 


cos 
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t 

Vi 

+  tz 

fi  = 
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sin  0 

7T 

"  26 

Then  the  integral  B  is  given  by  the  equation 


b 


a    C  +  t  1  -  (1  +  6)e~i 
6  J_b    "        ~5' 


B  =  nT  —^ wdt  (38) 


for  integration  along  both  horizontal  sides. 

Wherever  f <3 1  ^  1  the  integrand  is  replaced  by  the  infinite  series  in  the  equation 

1  -(1  +6)e~6       ~    .  ,    (-6)"  ,      , 

^-=t?o(*+1)(^  (39) 

Summation  of  this  series  is  continued  until  there  is  no  change  in  the  sum. 

GAUSS  INTEGRATION 

The  integration  is  completed  with  the  aid  of  16-point  Gaussian  integration 
multipliers.  Each  interval  of  integration  is  subdivided  into  subintervals  and  the  Gaussian 
integration  is  applied  progressively  to  successive  subintervals.  The  numbers  of 
subintervals  for  the  integration  of  A  and  B  are  the  integral  parts  of  the  expressions 

2+^J  2+^  (40) 

76  7a  V 

Then  the  Gaussian  integration  is  accurate  for  any  x  or  y.  For  large  values  of  x,  y  the 
time  for  computation  is  nearly  proportional  to  the  sum  of  |x|  and  \y\. 

Formulae  for  the  Taylor  series  expansion  of  the  source  pattern  at  nodal  points  in 
a  finite  grid  were  used  with  subroutine  RPLTM  to  obtain  gradients  at  several  grid  sizes. 
Extrapolation  to  infinite  grid  size  confirmed  the  line  integration  in  wave  number  space. 

The  components  of  the  gradient  of  potential  are  given  by  the  following  subroutine. 

SUBROUTINE  RPTNG   (AA,  AB,  AX,  AY,  AZ,   FX,   FY,   FZ) 

********************************************************************************** 

FORTRAN  SUBROUTINE  FOR  POTENTIAL  GRADIENT  OF  RECTANGULAR  PATTERN 

Httt4ttlittt*iHtt»Mt*tAt»(tU<lt(*tt(«(ttttttt»((llt(tltt*tt»*<ttttt*U*ttlt*tl****tttt*tt 

The  half-interval  a  of  the  pattern  is  given  in  argument  AA,  the  half-interval  6  of  the 
pattern  is  given  in  argument  AB,  and  the  Cartesian  coordinates  x,  y,  z  of  the  field  point 
are  given  in  arguments  AX,  AY,  AZ.  The  components  of  the  gradient  are  computed  by 
Gauss  integration  and  are  stored  in  the  functions  FX,  FY,  FZ. 


INTEGRATION  BY  PARTS 

For  the  integration  by  parts,  the  quadratic  formula  is  used  to  solve  the  equations 
for  6  to  give  t  as  a  function  of  6. 

For  the  vertical  sides  of  the  rectangle  the  variable  t  is  given  by  the  equation 


t 


xy{2-f  +  xx)  ±  |2|V(^  +  xx)2  -  (y2  +  z2) 


yz  +  zz 


Substitution  and  cancellation  lead  to  the  equation 


Vm*  = 


|z|(^  +  xx)  ±  iy^/(2-f  +  xx)2-(y2  +  z2) 


y2  +  z2 


and  differentiation  leads  to  the  equation 

xyV(2-f  +  xx)2-(y2  +  z2)  ±  \z\(^  +  ix) 


dt 

d~6 


2a 


(y2  +  z2)  L 


V(^  +  xx)2  -  (y2  +  e2) 


(41) 


(42) 


(43) 


The  ±  sign  in  the  equations  is  determined  by  substitution  of  the  limits  of  integration 
for  t.  The  sign  is  +  for  positive  t  and  is  -  for  negative  t. 

Substitution  in  the  equation  for  A  replaces  integration  with  respect  to  t  with 
integration  with  respect  to  <5.  The  natural  path  of  integration  in  the  £-plane  is  the 
real  axis,  whereas  the  natural  path  of  integration  in  the  6-plane  is  an  hyperbola  with 
vertex  where  <5  is  given  by  the  equation 


6  =  —  j|z| 
2a 


ix] 


(44) 


and  with  axis  in  the  direction  of  the  positive  real  axis.  There  is  a  singularity  on  each 
side  of  the  imaginary  axis.  Near  the  singularity  on  the  positive  side  6  is  given  by  the 
equation 


ca 


xx    +  e 


(45) 


Then  the  limiting  values  for  functions  at  the  singularity  are  given  by  the  equations 


xy 


V~vz  -  ~a 


VT+t2  = 


'y^  +  z' 


n/ 


yc  +  z' 


and  by  the  equation 


Ve 


dt 
d6 


<h/S 


(y2  +  z2)4 


(46) 


(e-0)  (47) 


The  hyperbolic  path  in  the  (5-plane  is  deformed  into  two  straight  lines  which  run  from 
the  positive  singularity  to  the  limits  of  integration  where  t  =  ±  a/6. 

For  the  horizontal  sides  of  the  rectangle  the  variable  t  is  given  by  the  equation 


t  = 


xx(2-¥  +  ^y)  ±  \z\^/(2-¥  +  xy)2-(x2  +  z2) 


(48) 


Substitution  and  cancellation  lead  to  the  equation 

v/rr7i  _  M(^  +  iy)  ±  ix%/(^  +  ry)*  -  (i»  4-  z') 


<49> 


and  differentiation  leads  to  the  equation 

(50) 


dt 


d6       (x2  +  Z2) 


izV^  +  xy)2  -  (x2  +  z2)  ±  \z\(*?  +  iy) 


V(^  +  iy)2-(x2  +  22) 


The  ±  sign  in  the  equations  is  determined  by  substitution  of  the  limits  of  integration 
for  t.  The  sign  is  +  for  positive  t  and  is  -  for  negative  t. 

Substitution  in  the  equation  for  B  replaces  integration  with  respect  to  t  with 
integration  with  respect  to  6.  The  natural  path  of  integration  in  the  i-plane  is  the 
real  axis,  whereas  the  natural  path  of  integration  in  the  6-plane  is  an  hyperbola  with 
vertex  where  6  is  given  by  the  equation 

6  =  ^{\z\-iy}  (51) 

and  with  axis  in  the  direction  of  the  positive  real  axis.  There  is  a  singularity  on  each 
side  of  the  imaginary  axis.  Near  the  singularity  on  the  positive  side  6  is  given  by  the 
equation 

6  =  ~  {Vx^W2  -  iy\  +  e  (52) 

do 

Then  the  limiting  values  for  functions  at  the  singularity  are  given  by  the  equations 

ix  i ~  \z\ 

t  =      ,  VI  +  i1  =      .  (53) 


Vx2  +  z2  Vx2 


+  2 


2 


and  by  the  equation 


Ve  —  -       '   '     "   3  (£  -  0)  (54) 

d6        (x2  +  22)4 


The  hyperbolic  path  in  the  (5-plane  is  deformed  into  two  straight  lines  which  run  from 
the  positive  singularity  to  the  limits  of  integration  where  t  =  ±  b/a. 
The  weighted  derivative 

is  expressed  as  a  power  series  expansion  by  11 -point  Lagrange  interpolation.  Let  6 
be  given  by  the  equation 

c5  =  77  +  e  (56) 

where  77  is  a  center  of  expansion  and  e  is  the  argument  of  expansion.  Inasmuch  as 
the  paths  of  integration  are  straight  lines  in  the  e-plane,  all  values  of  e  in  the  range 
of  interpolation  are  given  by  the  equation 

e  =U€U  (-llu^+1)  (57) 

where  eM  is  the  upper  limit  of  interpolation  and  u  is  the  variable  of  interpolation.  The 
limit  eu  can  be   cancelled  out  of  the  rational  expression  for  the  complex  Lagrange 
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coefficient  and  the  derivative  can  be  expanded  as  a  power  series  in  u.  In  a  preliminary 
computation  discrete  values  for  u  with  Chebyshev  spacing  and  discrete  coefficients 
were  computed  for  use  in  a  subroutine.  Thus  the  derivatives  are  expressed  as  power 
polynomials  in  the  ratio 

u=—  (58) 

while  division  by  eu  is  deferred  until  the  integration  by  parts. 

There  are  three  stages  of  expansion.  In  the  first  stage,  77  coincides  with  the  singularity 
and  derivatives  are  expanded  as  power  polynomials  in.e1/2.  The  range  of  interpolation 
is  extended  into  a  negative  range  in  order  to  obtain  data  for  central  interpolation. 
The  extension  in  the  e1/2-plane  corresponds  to  a  fold  in  the  e-plane.  The  interval  of 
integration  is  the  upper  half  of  the  range  of  interpolation.  In  the  later  stages  of 
expansion,  77  is  moved  further  out  along  the  path  of  integration  and  the  derivatives 
are  expanded  as  power  polynomials  in  e.  The  interval  of  integration  is  the  full  range 
of  interpolation. 

There  are  tight  limitations  on  the  range  of  interpolation.  The  radius  of  expansion 
near  the  first  singularity  must  be  a  small  fraction  of  the  distance  from  the  first 
singularity  to  the  second  singularity.  The  radius  of  expansion  near  the  second  center 
of  expansion  must  be  a  small  fraction  of  the  distance  from  the  center  of  expansion 
to  either  singularity. 

The  integration  of  each  series  in  powers  of  e1/2  is  completed  with  a  descending 
recurrence.  Required  for  the  integration  are  integrals  which  are  generated  by  the 
recurrence 


f 

Jo 


em  1 

eTn-ie-6de  =  — e~6  +  —   I    e^e'6  de  (59) 

771  77 


The  iteration  is  started  with  an  initial  approximation  for  771  =  32.  Completion  of 
integration  is  given  by  integrals  which  are  obtained  by  the  recurrence 

Jo  6  mr?  <5  m7?  Jo 

€m -7T2 de  (6°) 

where  771  is  alternately  half  an  odd  integer  and  half  an  even  integer.  The  iteration  is 
started  with  an  initial  approximation  for  771  =  32.  Stability  of  the  recurrence  is  achieved 
by  a  limitation  of  |e|  to  -7=  1*7 1- 

The  integration  of  each  series  in  powers  of  e  is  cycled  in  either  direction.  Required 
for  the  integration  are  integrals  which  are  generated  by  the  ascending  recurrence 

eme-6  de  =  -  eme~d  +  m       em-1e"*de  (61) 

Jo  t/0 

The  iteration  is  started  with  the  special  integral  in  the  equation 

e  e~6  de  =  e'*  -  (1  +  e)e~6  (62) 


J 

Jo 


and  the  iteration  is  continued  until  771  =  64.  If  |e|  ^  17  then  the  iteration  is  cycled  in 
descending  order  with  an  initial  approximation.  Completion  of  integration  is  given  by 
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integrals  which  are  generated  by  the  ascending  recurrence 

l-(l+<5)e-6  em      l-(l+<5)e-a  1 

Em — =— de  = I     eme_<5  de 

6  m  -  1  o  m  -  1 


77177 

771  -   1 


ff    m   .   1  -(1  +(5)e-6 


(63) 


where  m  is  a  sequence  of  integers.  The  iteration  is  started  with  the  aid  of  the  special 
integrals  in  the  equations 

l-(l+u)e-"    .         .       1-e-'  ~    (-6)'" 

„ 5? *"^ —g— £.W?sji  (64) 

6  1  -  (1  +  u)e"u  °°      (-6)k+z 

>_l_-L_*,.-i  +  .-.  +  7  +  toti-^-4,.EX-L_         (65) 

j  1  -  (1  +  u).-|d»  =  -  a  +  5  +  (2  +  d)e-'  —  £  (k  +("g78)t,  <^) 


I' 


and  the  iteration  is  continued  for  11  cycles.  If  |e|  ^  II77I  the  iteration  is  cycled  in 
descending  order  with  an  initial  approximation. 

In  the  actual  program  e  is  replaced  by  e/ew,  which  gives  integrals  ready  for  use 
with  the  Lagrange  coefficients,  and  scales  the  integrals  so  that  their  exponents  do  not 
exceed  the  tight  limitation  on  the  exponents  in  IBM  computers.  The  recurrence  equations 
can  be  verified  easily  by  differentiation 

Although  the  integration  by  parts  is  accurate  everywhere  that  \z\  is  not  small,  the 
integration  by  parts  breaks  down  when  yz  +  zz  is  zero  along  the  vertical  sides  and 
breaks  down  when  xz  +  zz  is  zero  along  the  horizontal  sides.  Under  such  circumstances 
the  program  falls  back  upon  Gauss  integration. 

For  large  values  of  x,  y  the  time  for  computation  is  independent  of  |x|  or  \y\.  The 
components  of  the  gradient  of  the  potential  are  given  by  the  following  subroutine. 

SUBROUTINE   RPTNX   (AA,  AB,  AX,  AY,  AZ,   FX,   FY,   FZ) 

******************************************************************************************* 

FORTRAN  SUBROUTINE  FOR  POTENTIAL  GRADIENT  OF  RECTANGULAR  PATTERN 

********************************************************************************************* 

The  half-interval  a  of  the  pattern  is  given  in  argument  AA,  the  half-interval  6  of  the 
pattern  is  given  in  argument  AB,  and  the  Cartesian  coordinates  x,  y,  z  of  the  field  point 
are  given  in  arguments  AX,  AY,  AZ.  The  components  of  the  gradient  are  computed  with 
integration  by  parts  and  are  stored  in  functions  FX,  FY,  FZ. 

GRADIENT  OVER  GRID 

Exploratory  computations  have  given  the  gradient  of  potential  on  a  vertical  line 
through  each  grid  point  of  the  source  pattern.  The  gradient  has  a  direction  which 
passes  through  the  center  and  has  a  magnitude  which  varies  with  height.  Far  from 
the  plane  the  magnitude  approaches  the  full  value  for  the  inverse  square  law,  while 
near  the  plane  the  magnitude  is  diminished  by  shielding  by  the  local  structure  of  the 
source  pattern.  As  distance  diminishes  toward  the  center  the  inverse  square  of  distance 
approaches  °°,  while  the  gradient  approaches  2tt. 
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The  potential  is  given  by  the  equation 

2ab    f+is  f+£  e-^Vl^«"^v) 
<p(x,y.z)  =  —   \  —       ,   „      _o       -da  d^ 


2b 


2a 


Va*  +  /?2 


(67) 


Combinations  of  the  partial  derivatives  replace  the  integrand  with  an  exact  differential 
which  can  be  integrated  at  once.  Results  of  integration  are  given  by  the  equation 


\       dx  dz )  TT      J  _JL 

2b 


2a 


dp 


(68) 


2a 


The  integral  is  zero  if  x  is  a  multiple  of  2a.  Results  of  integration  are  given  by  the 
equation 


dtp  d<p 

I  lz| y  — 

dy  dz 


2ab 


2a 


2a 


-N/a2+02|z|  +  t(ax+/Jy) 


2b 


da 


(69) 


2b 


The  integral  is  zero  if  y  is  a  multiple  of  26.  Thus  the  potential  gradient  has  a  radial 
direction  over  grid  points  where  x,  y  are  multiples  of  2a,  26.  The  integration  of  the 
inverse  square  over  the  surface  of  the  pattern  is  just  the  gradient  of  a  point  source 
when  the  field  point  is  over  a  grid  point. 

RECTILINEAR  LINE 

When  the  plane  contracts  to  a  line,  the  potential  and  the  potential  gradient  become 
symmetric  with  respect  to  the  line.  The  position  of  a  field  point  is  located  by  the 
cylindrical  polar  coordinates  r,  0,  z,  where  r  is  the  distance  from  the  line,  0  is  the 
azimuth  angle,  and  z  is  the  zenith  distance  along  the  line.  The  Fourier  transform  for 
a  function  on  the  line  is  given  by  the  equations 


A(k)  =  — 

V  27T 


f(z) 


f(z)e- 


dz 


A(K)e*"dK 


(70) 


(71) 


where  z  is  the  coordinate  in  physical  space  and  k  is  the  coordinate  in  wave  number 
space. 

A  universal  source  pattern  for  the  infinite  line  is  defined  by  the  equation 


sin 


/(*)  = 


7T2 

2c 


7T2 

2c 


(72) 


where  the  spacing  between  nodal  points  is  2c.  Application  of  the  Fourier  transform 
and  use  of  the  Euler  theorem  introduces  the  difference  of  two  cosines  and  the  sum 
of  two  sines  into  the  integrand  of  the  Fourier  amplitude.  The  difference  of  cosines  is 
an  even  function  of  z,  and  the  sum  of  sines  is  an  odd  function  of  z.  Both  functions 
are  divided  by  the  odd  function  z.  The  cosines  are  cancelled  by  the  integration.  There 
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remains  the  sum  of  sines  as  expressed  by  the  equation 

1     f+0°  sin(§  -  k)z  +  sin(£  +  k)z 

A(k)  =  — di 

4tt  J_m  lnz\ 


(73) 


2c 
Then  the  sines  cancel  each  other  whenever  /c  is  outside  the  range 


7T  7T  , 

iJSKS  +  ^  (74) 


and  the  amplitude  is  given  by  the  equation 

A(k)  =  -  (75) 

7T 

whenever  «;  is  inside  the  range.  Thus  the  function  f(z)  is  given  by  the  equation 

c    f  +  £     .  2c    f£ 

/(z)  =  -  etX2d/c  =  —  cos(kz)  d/c  (76) 

*J-£  ^    Jo 

where  the  range  of  integration  in  wave  number  space  is  finite. 

The  Laplacian  in  cylindrical  polar  coordinates  is  given  by  the  equation 

d*(p       1   dip        1    dz<p       d2<p  ,      x 

drc       r   dr       rc  90         dz 

The  solutions  of  Laplace's  equation  are  given  in  many  texts.  In  this  case  the  appropriate 
solution  contains  the  modified  Bessel  function  of  the  second  kind,  because  this  function 
diminishes  to  zero  with  increasing  argument2,3 
The  potential  is  given  by  the  equation 

2c    f  +  £     .  Ac    f£ 

tp(r,z)  =  —  etKZ  K0(kt)  die  =  —  cos(kz)  K0(kt)  die  (78) 

n   J-JL  tt     L 


where  K0(icr)   is   the   Bessel   function   of  order   zero.   The   Bessel   functions   satisfy  the 
equation 

—  K0(kt)  =  -kK^kt)  (79) 

dr 

where  Kl(icr)  is  the  Bessel  function  of  order  one.  The  Bessel  functions  satisfy  the  limits 
K0{Kr)  -*  -  log(/cr)  Kx{kt)  ->  —  (r  -  0)  (80) 

KT 

at  small  radius.  Thus  the  potential  satisfies  the  limit  equation 

-  27ir  -^  =  An  f(z)  (r-0)   (81) 

or 

where  f(z)  is  the  linear  density  a(z)  in  accordance  with  the  Gauss  theorem. 

13 


The  components  of  the  gradient  of  potential  are  given  by  the  equations 

It  IT 

=  +  —  k  e^'K^Kr)  dK  =  +  —  tc  cos(Kz)Ki(Kr)  die  (82) 

dr  7T   J_jl  rt   J0 


lie    f+zi        ■         /      x  4c    f^ 

k  etKZK0{Kr)  dK  =  +■ —  k  sin(icz)K0(icr)  dK 

7T      J_J1  7T     J0 


dip  2ic 

—  = I         k  eXK*K0{Kr)  dK  =  +  ■ —   I       k  sin(/c2)/f0(/cr)  dK  (83) 

dz  n    J_jl 

2c 


Inasmuch  as  these  integrals  are  limited  to  a  finite  range,  they  are  evaluated  with 
Gauss  integration  multipliers. 

The   potential  gradient   of   a   rectilinear   source   pattern   is  given   by   the   following 
subroutines. 

SUBROUTINE  BSSLK   (MO,  AZ,   IN,   FK) 
************************************************************************** 

FORTRAN  SUBROUTINE  FOR  MODIFIED  BESSEL  FUNCTION  OF  INTEGRAL  ORDER 

********************************************************************************************* 

The  mode  of  operation  is  given  in  MO.  The  real  and  imaginary  parts  of  the  argument  z 
are  given  in  array  AZ,  and  the  integer  order  n  is  given  in  IN.  The  modified  Bessel 
function  of  the  second  kind  is  computed  by  series  expansion,  rational  approximation, 
and  recurrence  relations.  If  MO  =  0,  the  real  and  imaginary  parts  of  the  function  Kn(z) 
are  stored  in  array  FK.  If  MO  =  1,  the  real  and  imaginary  parts  of  the  function  e*Kn(z) 
are  stored  in  array  FK. 

SUBROUTINE   LPTNG  (AC,  AR,  AZ,   FR,   FZ) 
********************************************************************************************* 

FORTRAN  SUBROUTINE  FOR  POTENTIAL  GRADIENT  OF  LINEAR  PATTERN 

********************************************************************************************* 

The  half-interval  c  of  the  pattern  is  given  in  argument  AC,  and  the  cylindrical 
coordinates  r,  z  of  the  field  point  are  given  in  arguments  AR,  AZ.  The  components  of 
the  gradient  are  computed  by  Gauss  integration  and  are  stored  in  the  functions  FR,  FZ. 

DISCUSSION 

The  ubiquitous  function 

sin  x 


(84) 
x 

is  sometimes  called  the  diffraction  function  because  it  appears  in  the  diffraction  of 
waves  at  an  aperture,  but  it  appears  in  other  applications  as  well.  It  should  have  a 
name  of  its  own.  It  is  proposed  that  it  be  called  the  sine  quotient  function. 

The  source  density  is  a  scalar  while  the  vortex  density  is  a  vector.  The  sine  quotient 
function  can  be  used  for  the  interpolation  of  the  source  density  or  either  of  the 
components  of  vortex  density.  Then  the  potential  gradient  of  the  pattern  function  is 
used  directly  with  source  density  in  the  computation  of  radial  flow,  or  in  a  vector 
product  with  either  component  of  vortex  density  in  the  computation  of  circular  flow. 

A  technique  for  the  computation  of  ideal  flow  over  arbitrary  bodies  has  been  invented 
by  Hess  and  Smith4,5.  The  same  technique  has  been  adopted  for  computation  of  flow 
over  ships  by  Dawson  and  Dean6,7.  In  this  technique  the  surface  of  the  body  is  subdivided 
by  a  grid,  then  each  panel  of  the  grid  is  projected  onto  a  plane  which  spans  the  four 
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corners  of  the  panel.  It  is  assumed  that  the  quadrilateral  projection  has  a  uniform 
source  density.  Each  plane  quadrilateral  is  assembled  from  plane  triangles,  for  which 
the  potential  gradient  is  given  by  analytical  expressions. 

In  the  case  of  an  infinite  plane,  uniform  quadrilaterals  would  be  reduced  to 
rectangular  plates.  The  potential  gradient  of  a  square  plate  is  compared  with  the 
potential  gradient  of  a  source  pattern  in  the  following  table. 

Potential  Gradients 
a  =  b  =  0.5  x  =  y  =  0 


2 

Inverse 
Square 

Plate 

Pattern 

0 

oo 

2tt 

2tt 

0.0625 

256 

5.58063622 

5.41516931 

0.125 

64 

4.90438059 

4.68192879 

0.25 

16 

3.70918087 

3.53400628 

0.5 

4 

2.09439510 

2.09491828 

1.0 

1 

0.80543168 

0.86196571 

2.0 

0.25 

0.23543002 

0.24793473 

3.0 

0.11111111 

0.10812127 

0  11106356 

4.0 

0.0625 

0.06154089 

0.06249869 

5.0 

0.04 

0.03960461 

0  03999996 

6.0 

0.02777777 

0.02758643 

0  02777777 

7.0 

0.02040816 

0.02030466 

0  02040816 

8.0 

0.015625 

0.01556424 

0  01562500 

The  entries  in  the  table  were  computed  with  the  aid  of  RPLTG,  RPTNG,  and  RPTNX.  The 
table  shows  that  the  inverse  square  law  gives  accurate  gradients  for  z  >  4. 

For  rectangular  plates  the  field  point  must  be  kept  away  from  the  edges,  where 
the  gradient  is  infinite,  whereas  for  pattern  functions  there  is  no  restriction  on  the 
location  of  the  field  point. 

CONCLUSION 

Pattern  functions  which  are  based  upon  the  sine  quotient  function  are  useful  for 
the  computation  of  the  potential  and  the  potential  gradient  of  source  distributions 
in  the  plane.  Directly  over  a  grid  point  the  potential  gradient  has  a  direction  which 
passes  through  the  center  of  the  pattern,  but  has  a  magnitude  which  is  shielded  from 
the  full  inverse  square  law. 
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